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Abstract 


This  report  presents  the  results  of  investigations  in  using  a  medium-to-high  accuracy  inertial 
navigation  system  (INS)  and  a  geodetic  quality  Global  Positioning  System  (GPS)  receiver  on  an 
airborne  platform  to  determine  the  gravity  disturbance  vector  along  the  trajectory  of  the  aircraft. 
The  main  products  of  these  investigations  consist  of  derivations  of  alternative  algorithms  to 
estimate  the  horizontal  gravity  vector  components  given  that  the  INS’s  under  consideration  have 
different  types  of  available  output  data,  either  time-integrated  accelerations  and  direction  cosines  or 
raw  accelerometer  and  gyro  data.  In  addition,  two  distinct  data  processing  methodologies  are 
elucidated  to  combine  the  GPS  and  INS  data.  The  first  is  based  on  a  traditional  Kalman  filter 
approach  that  treats  the  gravity  disturbance  vector  as  a  stochastic  process,  and  the  second  is  an 
extension  of  the  scalar  airborne  gravimetry  concept.  Also,  the  configuration  details  of  a  loosely 
coupled  INS/GPS  system  for  airborne  application  are  given,  with  particular  emphasis  on  time 
synchronization. 

A  limited  number  of  field  tests  (ground  and  airborne)  with  the  LN93  and  LNIOO  were 
performed.  Results  from  these  tests  verify  the  specified  accuracy  of  the  INS  and  yield  initial 
attempted  estimations  of  the  vertical  deflection  using  the  traditional  Kalman  filter  approach.  While 
the  latter  are  inconclusive  in  assessing  the  precision  level  achievable  with  this  method  and 
instrumentation  (due  to  unexpectedly  poor  gyro  performance),  they  do  indicate  the  feasibility  of  the 
approach  provided  the  INS  gyros  perform  at  the  specified  level. 


1  Background 


Knowledge  of  the  Earth's  gravity  field  is  an  essential  element  in  most  geodetic  and  many 
geophysical  applications,  dating  back  to  the  17th  century  (Newton,  Huygens),  and  firmly 
established  for  geodesy  by  F.R,  Helmert  (1880)  who  developed  the  modem  concept  of  the  geoid 
as  geodetic  reference  surface  to  approximate  the  shape  of  the  Earth.  (The  geoid  is  an  equipotential 
surface  of  the  Earth’s  gravity  field  that  closely  approximates  mean  sea  level.)  Since  the  1950’s,  the 
Earth’s  gravitational  field  also  played  a  significant  role  in  satellite  orbit  determination  and  in  inertial 
navigation  and  guidance.  Conversely,  with  observations  of  satellite  orbits,  and  especially  with 
satellite  altimetry  over  the  oceans  (approximately  the  geoid),  as  well  as  the  steady  accumulation  of 
surface  gravity  measurements,  very  accurate  models  of  the  Earth’s  gravity  field  have  been 
developed  during  the  last  forty  years,  and  significant  strides  have  been  made  during  the  last  decade 
(Nerem,  1995;  Nerem  et  al.,  1995).  With  these  efforts,  the  geoid  can  now  be  calculated 
worldwide  to  an  accuracy  of  better  than  one  meter  (NIMA,  1997);  in  ocean  areas  to  about  30  cm 
(Rapp,  1996),  and  on  some  continents  to  better  than  2  cm  over  50  km  (Milbert,  1997;  see  also 
Balmino  and  Sanso,  1995).  Yet  the  gravity  field  with  its  infinite  detail  and  information  remains  a 
key  component  in  geodesy  and  geophysics  and  the  object  of  continually  improving  measurement 

One  particularly  relevant  need  to  know  the  Earth’s  gravity  field  concerns  inertial  navigation  and 
guidance  systems.  Historically,  the  inertial  navigation  system  (INS)  has  provided  the  principal 
means  of  autonomous  navigation  and  guidance  systems  of  aircraft,  ships  (including  submarines), 
and  missiles.  The  recent  rapid  developments  and  achievement  of  high  accuracy  with  the  Global 
Positioning  System  (GPS)  have  quickly  overtaken  the  INS  in  terms  of  the  range  of  applications 
and  the  long-term  maintenance  of  mission  accuracy.  However,  due  to  the  recognition  of  the 
limitations  and  vulnerabilities  of  GPS,  as  well  as  new  technological  advancements  in  inertial 
instrumentation,  INS  is  considered  a  valuable,  if  not  indispensable  means  to  enhance,  aid,  or 
replace  GPS  in  various  navigation,  positioning,  and  guidance  scenarios.  Inertial  navigation  still  is 
the  only  method  that  provides  fully  autonomous,  completely  stealthy,  and  jam-proof  navigation 
and  guidance.  Efforts  to  make  GPS  immune  to  jamming  have  been  largely  unsuccessful. 

On  the  other  hand,  unlike  GPS,  very  precise  inertial  navigation  systems  require  knowledge  of 
the  Earth’s  local  gravity  vector  field,  specifically  the  direction  of  the  gravity  vector.  It  may  be 
noted  that  GPS  also  requires  some  knowledge  of  the  Earth  gravitational  field  in  order  to  determine 
precise  orbits  of  the  GPS  satellites;  however,  only  very  long-wavelength  information  is  required 
since  the  orbits  have  high  altitude  (20,000  km).  The  direction  of  gravity  is  prescribed  naturally  by 
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the  total  mass  density  structure  of  the  Earth  and  can  be  referenced  to  the  direction  deHned  by  the 
surface  normal  (perpendicular)  of  a  rotational  ellipsoid  that  approximates  the  geoid.  The  difference 
between  these  two  directions  is  known  as  the  deflection  of  the  vertical.  The  two  angular 
components  of  the  deflection  of  the  vertical  are  also  related  directly  to  the  horizontal  (with  respect 
to  the  ellipsoid)  components  of  the  gravity  vector  (Figure  1),  where  the  signs  are  a  matter  of 
convention. 


P 


Figure  1:  The  geometry  of  the  deflection  of  the  vertical. 

Although  global  models  for  the  Earth's  gravitational  field  are  improving  (like  NIMA’s  EGM96, 
Lemoine  (1996))  shorter-wavelength  signatures  are  either  very  poorly  modeled  Oess  than  100  km) 
or  only  moderately  well  known  (100  km  to  1000  km),  especially  in  rugged  terrain  regions  of  the 
world  and  the  polar  regions.  A  recent  study  by  Jekeli  (1998a),  that  compares  astronomically 
determined  deflections  of  the  vertical  in  the  United  States  to  corresponding  deflections  computed 
from  EGM96,  indicates  that  the  EGM96  gravity  model  at  wavelengths  shorter  than  200  km  may  be 
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under-powered.  This  conclusion  (requiring  further  verification)  is  indicated  from  considerations  of 
the  implied  power  of  the  signal  beyond  the  cutoff  wavelength  of  100  km.  Specifically,  the  EGM% 
model  implies  that  the  root-mean-square  of  the  deflection  of  the  vertical  with  wavelengths  shorter 
than  100  km  in  the  conterminous  U.S.  (CONUS)  is  only  1.33  arcsec;  while  the  corresponding 
value  for  the  astronomically  determined  deflections  in  CONUS  is  3.84  arcsec. 

On  the  other  hand,  a  significant  part  of  the  error  budget  of  high  accuracy  inertial  systems  (0.1 
to  1.0  nmi/hr  systems)  is  the  unknown  gravity  field,  as  shown  by  Jekeli  (1997).  That  is,  a  high 
accuracy  system  which  should  be  accurate  to  0. 1  nautical  mile  per  hour  (or,  assuming  linearity, 
about  15  meters  per  5  minutes)  is  susceptible  to  unknown  gravity  effects  (even  with  respect  to 
EGM96)  of  about  10  meters  over  that  same  interval.  Also,  it  is  seen  in  this  study  that  gravity 
compensation  for  high-accuracy  INS  is  required  primarily  for  wavelengths  greater  than  10  km,  or 
higher  depending  on  the  speed  of  the  vehicle.  Therefore,  better  knowledge  of  the  Earth’s  gravity 
field,  especially  at  wavelengths  between  10  km  and  200  km,  is  required  to  support  high  accuracy 
(free-)  inertial  worldwide  navigation  and  guidance. 

Lacking  accurate  and  sufficiently  dense  (vertical)  gravity  surveys,  the  total  gravity  vector  is  still 
not  well  known  in  remote  and  inaccessible  regions,  such  as  the  Himalayas,  the  Andes,  the  polar 
regions,  some  coastal  regions,  and  some  interior  continental  regions.  To  support  improved 
modeling  of  the  gravity  vector  field  requires  regional  gravity  mapping.  The  more  direct  method 
using  astronomic  measurements  to  determine  the  direction  of  gravity  is  considered  too  costly,  as 
well  as  practically  unfeasible,  in  many  regions.  However,  it  is  possible  to  estimate,  on  the  basis  of 
potential  theory,  the  horizontal  gravity  components  from  vertical  component  values  that  are  densely 
distributed  regionally.  Therefore,  to  date,  most  determinations  of  the  deflections  have  been  derived 
from  scalar  (vertical)  gravimetry.  That  is,  only  the  vertical  component  of  the  gravity  vector  is 
measured  since  this  is  much  easier  in  terms  of  methodology  and  cost.  These  surveys  can  be  done 
on  the  ground  using  conventional  gravimetric  surveys  or  with  much  greater  efficiency  and 
accessibility  to  remote  areas,  but  with  less  accuracy,  using  airborne  platforms.  For  example,  one 
of  the  most  successful  of  the  airborne  gravimetric  systems  is  being  operated  by  the  Naval  Research 
Laboratory  (NRL)  by  Dr.  John  Brozena  (Brozena  and  Peters,  1994).  Recently,  vast  previously 
unmq)ped  areas  of  Greenland  and  the  Arctic  Ocean  were  overflown  with  a  gridded  pattern  by  NRL 
yielding  an  accuracy  of  about  5  mgal  with  15  km  resolution  in  the  vertical  gravity  values  (Forsberg 
and  Kenyon,  1994). 

A  significant  drawback  of  scalar  surveys  is  that  the  determination  of  the  horizontal  gravity 
components  requires  numerical  integration  of  gravity  data  over  larger  regions  (Vening-Meinesz 
integrals).  Thus,  mapping  the  total  gravity  vector  using  scalar  gravimetry  relies  not  only  on 
accurate  measurements  over  large  regions,  but  also  on  a  physical  model  (usually  approximate)  that 
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relates  the  measurements  to  the  desired  quantities,  in  this  case,  horizontal  gravity  components. 

Other  moving-base  systems  that  measure  the  total  gravity  vector  or  the  gravity  gradient  tensor 
have  been  considered  (and  tested,  in  the  case  of  the  gradiometer,  Jekeli  (1988)).  These  systems 
provide  a  more  direct  (in  situ)  measure  of  the  vector  field  and  rely  less  on  the  use  of  approximate 
potential  theory  models  to  estimate  the  deflection  components.  A  fully  operational  system,  the 
Gravity  Gradiometer  Survey  System  (GGSS)  constructed  by  Bell  Aerospace,  was  tested  under 
NIMA  (formerly  DMA)  sponsorship  in  the  1980’s  and  yielded  3-6  mgal  accuracy  of  wavelengths 
between  10  and  1(X)  km  in  the  horizontal  components  of  the  gravity  vector  (Jekeli,  1993).  This 
particular  gravity  gradiometer  system  measures  the  five  independent  gravity  gradients  that  can  be 
integrated  to  yield  the  gravity  vector.  The  gradiometer  is  considered  in  the  geophysical  community 
to  be  the  instrument  that  will  yield  the  high  resolution  and  accuracy  (1  mgal  over  1  km)  gravity 
surveys  needed  for  mineral  exploration,  upper  crustal  structure  modeling,  and  local  geoid 
determination.  Therefore,  gravity  gradiometry  promises  to  yield  directly  the  fine  structure  of  the 
vector  gravity  field,  while  also  overcoming  the  scalar  limitations.  Further  development  of 
instruments  and  testing  are  needed  to  demonstrate  the  full  potential  of  these  systems. 

A  third  system  that  has  been  under  study  for  some  time,  but  not  fully  realized  in  practice,  is  the 
integrated  INS/GPS  platform.  The  principle  of  measurement  is  identical  to  moving-base  scalar 
gravimetry,  but  extended  to  three  dimensions.  With  GPS  one  measures  the  total  kinematic  three- 
dimensional  acceleration  of  the  vehicle  as  determined  from  its  three-dimensional  positions.  Using 
the  triad  of  accelerometers  in  the  INS  one  measures  the  three-dimensional  specific  forces  acting  on 
the  vehicle.  The  difference  between  the  two  types  of  acceleration  vectors  is  the  gravitational 
vector.  As  such,  the  INS/GPS  system  provides  direct  measurements  of  the  horizontal  (and 
vertical)  gravity  components  that  rely  only  on  the  quality  of  the  measurement  systems  and  not  on 
(usually  approximate)  analytic  models  relating  vertical  and  horizontal  components.  On  the  other 
hand,  the  quality  of  the  sub-systems,  INS  and  GPS,  is  such  that  errors  at  both  ends  of  the 
spectrum  limit  the  gravitational  signal  that  may  be  extracted  from  such  an  integration. 
Nevertheless,  numerous  studies  have  been  conducted  that  show  the  feasibility  of  integrating  INS 
with  GPS  to  conduct  airborne  vector  gravimetry  (Eissfeller  and  Spietz,  1989;  Schwarz  et  al.,  1992; 
Jekeli,  1995).  Simulation  studies  by  Jekeli  (1995),  Gleason  (1992),  and  Wang  et  al.  (1997), 
among  others,  have  shown  that  the  accuracy  of  recoverable  gravity  is  about  3-4  mgal  over 
wavelengths  of  30  km  to  200  km  using  a  high-accuracy  INS  and  differential  GPS  positioning. 
The  development  of  such  systems  has  been  hindered  in  part  by  the  expense  of  high  quality  INS, 
being  about  five  times  the  cost  of  a  high  quality  GPS  receiver. 

The  most  difficult  aspect  of  vector  gravimetry  using  GPS/INS  is  the  drift  error  in  the  angular 
data  provided  by  the  gyros  that  are  needed  to  properly  orient  the  accelerometers.  An  orientation 
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error  at  the  level  of  1  arcsec  translates  into  a  5  mgal  error  in  the  horizontal  gravity  component  The 
growth  in  the  orientation  error  can  be  controlled  by  imposing  velocity  control  (zero  velocity 
updates,  or  ZUPT’s)  at  frequent  intervals  (e.g..  Rose  and  Nash,  1972);  or  on  an  airborne  system 
by  using  a  star  tracker  to  provide  absolute  orientation  of  the  system  to  an  external  reference  frame, 
as  was  done  with  success  by  Northrop  Corporation  (1986).  Alternative,  less  costly  schemes  have 
been  considered,  such  as  photogrammetric  control  (Wang  et  al,  1997).  Or,  Jekeli  (1995)  also 
argued  that  if  the  orientation  error  growth  is  indeed  a  long  wavelength  phenomenon,  then  short- 
wavelength  gravity  estimation  should  be  possible  without  resorting  to  external  attitude  aiding, 
where  the  long-wavelength  component  (>  200  km)  would  come  from  a  global  spherical  harmonic 
model. 

This  report  reviews  some  of  the  activities  surrounding  the  use  of  the  LN93  and  LNIOO  Litton 
inertial  navigation  systems  and  the  Global  Positioning  System  (GPS)  to  measure  the  total  gravity 
vector  using  an  airborne  platform.  No  conclusive  results  from  actual  test  flight  can  be  reported  at 
this  time;  however,  various  models  and  evaluations  of  the  INS  are  presented,  with 
recommendations  for  further  studies. 


2  INS  Instrumentation  and  Measurements 


An  inertial  navigation  system  (INS)  consists  of  an  inertial  measurement  unit  (IMU)  and  a 
navigation  computer.  The  essential  element  of  an  IMU  is  the  accelerometer  whose  output  is 
integrated  twice  in  time  to  obtain  positions.  A  common  class  of  accelerometers  is  the  force- 
rebalance  type.  Although  a  variety  of  designs  exist,  most  are  based  on  the  principle  of  maintaining 
the  null  position  of  a  proof  mass  on  a  spring.  The  electronicaUy  applied  force  needed  to  do  this  is  a 
measure  of  the  acceleration.  Three  accelerometers  with  sensitive  axes  mutually  perpendicular 
provide  three-dimensional  navigation. 

Of  equal  importance,  however,  is  the  coordinate  frame  in  which  the  accelerometers  are  to 
provide  positions.  In  this  respect,  not  only  the  orientation  of  the  accelerometers  in  the  coordinate 
system,  but  also  their  angular  velocity  affect  the  determination  of  position.  The  orientation 
determines  the  component  of  the  position  vector  that  a  particular  accelerometer  provides;  and,  as  is 
known  from  elementary  physics,  the  angular  rates  with  respect  to  an  inertial  frame  contaminate  the 
accelerometer  outputs  with  centrifugal  and  Coriolis  accelerations. 

The  orientation  and  angular  rates  of  the  accelerometer  platform  are  determined  with 
gyroscopes.  Most  INS  for  commercial  deployment,  today,  use  either  the  ring  laser  gyro  or  the 
less  costly  (and  less  accurate)  fiber  optic  gyro.  Both  types  of  gyros  require  that  the  INS  be 
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mechanized  in  the  so-called  strapdown  configuration;  that  is  the  INS  is  physically  mounted  to  the 
frame  of  the  vehicle  (no  gimbal  support  system).  This  reduces  the  cost  of  the  INS  considerably  as 
compared  to  gimbal-supported,  local-level  stabilized  systems.  It  also  subjects  the  system  to  the 
entire  spectrum  of  vehicle  dynamics  and  the  navigation  accuracy  is  less  than  with  the  stabilized 
(motion-isolated)  systems.  Figure  2  is  a  schematic  drawing  of  a  strapdown  platform  for  the 
IMU’s. 


The  ring  laser  gyro  (RLG),  in  principle,  has  no  moving  parts  -  it  is  based  on  the  fact  that  the 
frequency  of  a  light  beam  travelling  in  a  resonant  closed  circuitous  path  (where  the  number  of 
wavelengths  is  always  the  same;  i.e.,  it  is  a  laser)  must  change  if  the  apparent  length  of  the  circuit 
changes  because  it  rotates  in  inertial  space.  The  latter  property  is  called  the  Sagnac  effect  and  is  a 
consequence  of  general  relativity  where  the  speed  of  light  is  a  constant  independent  of  the  frame  in 
which  it  is  traveling.  Two  such  counter-travelling  beams  of  light  are  used  to  create  a  fringe  pattern 
where  they  recombine.  The  fringe  pattern  is  stationary  if  there  is  no  rotation  in  inertial  space. 
But,  it  migrates  in  the  presence  of  rotation  about  the  axis  perpendicular  to  the  plane  of  the  circuit, 
because  then  one  beam  sees  a  longer  path,  the  other  a  shorter  path;  and  the  number  of  fringes 
passing  a  detector  per  unit  time  indicates  the  rate  of  rotation. 

The  major  problem  with  RLG's  is  a  phenomenon  called  lock-in:  due  to  imperfections 
(scattering  of  light  in  the  resonant  cavity)  the  two  beams  lase  at  the  same  frequency  even  in  the 
presence  of  a  small  rotation,  typically  up  to  several  hundred  degrees  per  hour.  They  lock  to  the 
same  frequency  and  indicate  a  zero  rotation.  One  common  solution  to  this  problem  is  to  bias  the 
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output  of  the  device  by  applying  a  physical  rotation  away  from  the  lock-in  range.  To  maintain 
stability,  this  mechanical  "bias"  is  in  the  form  of  an  alternating  rotation,  i.e.,  a  dithering  or 
oscillation  (tens  to  hundreds  of  Hz)  of  the  gyro  about  its  sensitive  axis.  This  is  the  operational 
concept  of  the  gyros  used  in  the  Litton  LN93  INS. 

A  newer  approach  is  to  apply  an  optical  bias  by  creating  left-  and  right-circular  polarization  of 
the  two  beams,  respectively.  In  the  presence  of  an  applied  magnetic  field,  the  speeds  of  the  beams 
differ,  which  is  equivalent  to  an  effective  difference  in  path  length  (Faraday  effect),  hence  creating 
a  bias  in  the  frequency  difference.  The  complete  absence  of  moving  parts  (no  mechanical 
dithering)  improves  the  stability  and  substantially  reduces  the  random  noise  of  the  IMU  as  a  whole. 
This  is  the  operational  concept  of  the  gyros  used  in  the  Litton  LNIOO  INS. 

Size,  weight,  and  power  requirements,  as  well  as  the  standard  error  budgets  for  each  system 
are  listed  in  Table  1  and  Table  2.  The  accelerometers  are  the  same  model  for  both  systems,  but 
presently  manufactured  sensors  are  slightly  improved.  The  characteristics  of  these  two  systems 
were  provided  by  Litton  (person  communication).  Both  are  quite  similar  and  several  other 
subordinate  error  sources,  especially  misalignments,  temperature  transients,  and  other  correlated 
noise,  as  well  as  acceleration  sensitivities  are  not  included  here. 


Table  1:  Essential  characteristics  of  the  LN93  and  LNIOO  inertial  systems. 


LN93 

LNIOO 

Data  Rate 

20  Hz  (user  defined) 

32  Hz  (256  Hz,  raw  data) 

Size  (excludes  mount) 

1089  cu.  in. 

539  cu.  in. 

Weight  (excludes  mount) 

48.5  lbs. 

19.4  lbs 

Power 

28VDC,  150  Watts 

28  VDC,  26.5  Watts 
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Table  2:  Essential  error  budget  for  LN93  and  LNIOO  inertial  systems. 


LN93 

LNIOO 

Accelerometer  (Litton  A-4  model) 

Bias  Error 

25mGal 

20  mGal 

Scale  Factor  Error 

120  ppm 

40  ppm 

White  Noise 

5  mGal/VHz 

5  mGal/VHz 

Gyro  (ring  laser  gyro) 

Bias  (Drift)  Error 

0.003  °/hr 

0.003  °/hr 

Scale  Factor  Error 

5  ppm 

0.2  ppm 

White  Noise 

0.0015  7Vhr 

<0.001  7Vhr 

Figure  3  shows  the  basic  configuration  of  the  INS  and  computer  interface  for  the  LN93.  A 
similar  setup  exists  for  the  LNIOO.  The  computer  used  with  the  LN93  is  a  laptop  computer 
mounted  in  a  so-called  docking  station  that  contains  both  the  1553  bus  controller  card  and  the  timer 
card.  The  1553  card  is  the  direct  interface  to  the  INS  and  all  requests  of  data  from  the  INS  are 
made  through  this  bus  controller.  The  timer  card  can  be  programmed  to  synthesize  a  sequence  of 
interrupts  at  a  specified  rate  which  is  then  used  to  request  the  INS  data  at  that  rate.  In  the  case  of 
the  LN93,  the  data  request  rate  is  20  Hz.  All  INS  data  requests  are  synchronized  to  GPS  time. 
Each  data  item  from  the  INS  includes  a  time  tag  relative  to  the  INS  clock  which  can  be  used  to 
determine  the  actual  time  for  which  the  data  item  corresponds.  Additional  data  are  retrieved  from 
the  GPS  receiver  by  the  bus  controller  software  to  determined  the  actual  GPS  time  stamp  of  the 
one-second-pulse,  and  to  provide  altitude  data  for  the  INS.  The  INS  data  are  saved  by  the  bus 
controller  software  in  the  RAM  of  the  laptop  computer  and  transferred  to  hard  disk  storage  at  the 
conclusion  of  the  mission.  The  GPS  raw  data  are  stored  in  the  GPS  receiver  and  also  later 
retrieved  and  processed. 


Figure  3:  Basic  configuration  of  INS,  GPS  receiver,  and  computer  test  equipment 
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The  data  output  options  of  an  INS  range  from  position  and  velocity  to  orientation  angles,  as 
well  as  raw  sensor  data,  system  status  data,  and  time  tags.  The  type  of  output  is  entirely  dictated 
by  the  navigation  software  that  the  INS  manufacturer  creates  for  its  customers.  The  software  for 
the  LN93,  that  is  on  loan  to  Ohio  State  University  from  Litton,  allows  users  to  request  position  and 
velocity  in  the  navigation  frame  (wander-azimuth  frame)  and  the  navigation-frame-to-Earth-fixed- 
frame  direction  cosine  mauix,  all  with  double  precision  (32  bit  data),  while  orientation  angles  of 
the  body  frame  with  respect  to  the  navigation  frame  are  available  only  in  single  precision  (16  bit 
data).  No  raw  accelerometer  and  gyro  data  are  available.  The  position  coordinates  have 
quantization  errors  of  about  43  cm  in  latitude  and  63  cm  in  longitude,  being  due  entirely  to  an 
artifact  of  the  data  processing  performed  by  the  navigation  software.  This  is  shown  in  Figure  4 
where  differences  between  indicated  latitude  and  longitude  values  (converted  to  linear  measures  in 
units  of  meters)  where  computed  for  one-second  intervals.  The  data  were  obtained  by  moving  the 
LN93  in  the  laboratory.  The  velocities  do  not  have  this  quantization  error  and  can  be  integrated 
independent  of  the  INS  software  to  obtain  position  coordinates  at  a  precision  level  of  a  centimeter, 
or  better. 

Nevertheless,  the  lack  of  availability  of  raw  sensor  data  from  the  LN93  makes  it  less  useful  for 
wider  geodetic  applications  and  tightly  coupled  integration  with  GPS.  In  addition,  care  must  be 
exercised  in  synchronizing  the  INS  and  GPS  data,  since  the  INS  software,  performing  real-time 
data  processing,  makes  data  available  to  the  user  with  time  validity  tags  showing  significant  time 
delays  from  the  actual  data  output  time.  The  software  of  the  LNIOO  obtained  by  OSU’s  Center  for 
Mapping  was  specifically  modified  to  provide  raw  sensor  data  which  can  be  processed  by  the  user 
independently  of  the  navigation  software.  Furthermore,  the  time  synchronization  is  less 
problematic  since  the  time  of  raw  sensor  output  coincides  very  closely  to  the  time  validity  of  the 
data  (since  very  little  processing  takes  place  in  the  INS).  The  raw  data  for  the  LNIOO  consist  of 
velocity  increments  from  the  accelerometers  and  angle  increments  from  the  gyros  at  the  high 
sampling  rate  of  256  Hz.  In  Section  3,  algorithms  are  developed  to  process  these  data  for  geodetic 
applications. 
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3  Algorithms  and  Models  for  Vector  Gravimetry 

3. 1  Basic  Equations 

Navigation  using  an  INS  requires  knowledge  of  the  Earth’s  gravity  vector  field.  That  is,  the 
INS  operates  under  the  basic  relationship  (Newton’s  Second  Law  of  Motion)  between  specific 

force  sensed  by  accelerometers,  a ,  and  total  acceleration  of  the  vehicle,  the  second  time  derivative 
of  position,  x ,  given  by: 

x  =  a  +  g  (jj 

where  g  is  the  gravitation  vector.  Integrating  this  equation  to  get  position,  x,  and  velocity,  x, 
from  sensed  acceleration,  a,  requires  knowledge  of  the  gravitation  vector,  g.  Conversely, 
combining  INS  with  an  independent  measurement  of  x  (or  x ),  such  as  from  GPS,  yields  a 
measurement  of  g . 

Equation  (1)  holds  for  an  inertial  frame,  only.  For  frames  attached  to  the  rotating  Earth,  or  to  a 
vehicle,  similar  equations  are  readily  derived,  but  are  complicated  by  the  inclusion  of  Coriolis 
terms  arising  from  the  rotation  of  the  Earth  or  vehicle  with  respect  to  an  inertial  frame  (see  Britting, 
1971;  Jekeli,  1998b).  Usually,  for  local-level  mechanizations,  the  gravitation  vector  is  combined 
with  the  centrifugal  acceleration  arising  from  Earth’s  rotation  and  is  then  called  the  gravity  vector. 
For  the  sake  of  simplicity  we  use  the  notation,  g ,  for  the  gravity  vector,  as  well,  since  the 
gravitation  vector  will  not  reappear  in  our  discussions.  The  gravity  vector  may  be  separated  into 
two  parts,  a  dominant,  so-called  normal  gravity  vector,  y,  that  accounts  for  the  centrifugal  part 
and  all  but  about  one  part  in  lO’*  of  the  total  gravitational  vector,  and  a  residual,  the  so-called 
disturbance  gravity  vector,  6g ,  such  that  g  =  y  n  +  6g ,  where  n  is  a  unit  vector  along  the  normal 
to  an  ellipsoid  and  y  =  |  y|  .  The  value  y  can  be  calculated  exactly  on  the  basis  of  a  well  defined 
normal  gravity  field  (Heiskanen  and  Moritz,  1967),  where  for  highest  accuracy  at  points  above  the 
ellipsoid,  the  curvature  of  the  field’s  horizontal  gradient  must  be  taken  into  account.  In  a  local 
north-east-down  coordinate  system 

5g  =  g-Y=(g4  -gil  5g)^  (2) 

(see  Figure  1)  where  ^ ,  q  are  the  deflection  of  the  vertical  components,  and  6g  =  |  g  |  -  y  is  the 
(scalar)  gravity  disturbance. 

The  analogue  to  equation  (1),  now  formulated  for  the  north-east-down  coordinate  frame,  is 
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written  as 


^v"  =  a"-(cD? +Q)P)xv"  +  g"  (3) 

where  v"  is  the  velocity  in  the  navigation  frame  (n-frame);  a"  and  g"  are,  respectively,  the 
sensed  acceleration  and  gravity  vector,  both  coordinatized  in  the  n-frame;  and  ea?, ,  <0?  are 
angular  rate  vectors,  respectively,  of  the  n-frame  and  of  the  Earth-fixed  frame  (e-frame)  each  with 
respect  to  the  inertial  frame  (i-frame),  as  coordinatized  in  the  n-frame  (Britting,  1971).  The  second 
term  on  the  right  of  (3)  is  essentially  due  to  the  Coriolis  and  centrifugal  accelerations  of  the  n-frame 
with  respect  to  the  i-frame. 

The  LN93  navigation  software  uses  the  so-called  wander-azimuth  frame  to  perform  the 
integration  of  the  navigation  equations.  This  avoids  certain  singularities  in  the  navigation  equations 
near  the  polar  regions  where  the  meridians  converge  (and  longitude  rate  becomes  undefined).  The 
wander  azimuth  frame  is  related  to  the  n-frame  as  shown  in  Figure  5.  The  navigation  equations  in 
this  case  are  given  by 

iv»  =  a’*’-«  +  <)xv»  +  g»'  (4) 

where  v'*'  is  the  velocity  in  the  w-frame,  and  a'*' ,  g^  similarly  are  coordinatized  in  the  w-frame. 
The  sum  of  angular  rate  vectors  in  (4)  is  also  equal  to 

(5) 

=  (B^^  +  2C^CD? 


where  is  the  transformation  matrix  from  the  e-frame  to  the  w-frame.  It  is  determined  by 
integrating  the  angular  rates,  ,  given  approximately  by 


W 

ew 


f  v^ 

-vr 

0  ; 


(6) 


where  r  is  the  geocentric  radius  to  the  vehicle.  The  details  of  this  integration  are  analogous  to  the 
procedure  outlined  below  for  the  raw  gyro  data  processing. 
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Figures:  Wander-azimuth  frame 

As  a  first  approximation,  navigation  computers  use  only  the  readily  computable  normal  gravity, 
Y"  =  =  Y  n ,  in  (3)  or  (4).  The  velocity  indicated  by  the  LN93  INS  is  then  given  (erroneously) 

by  v* ,  which  satisfies  the  differential  equation: 

=  +  (7) 

That  is,  this  is  the  equation  that  is  integrated  by  the  LN93  navigation  computer  to  obtain  the  w- 

frame  velocities.  The  latitude  and  longitude  can  be  extracted  directly  from  the  transformation 
matrix,  ,  given  by 


pw  _ 


-  cosa  sin<f»  cosX  -  sina  sinX  -  cosa  sin<{)  sinX  +  sina  cosX  cosa  cos<|> 
sina  sin<|)  cosX  -  cosa  sinX  sina  sin<J>  sinX  +  cosa  cosX  -  sina  cos<{> 
cos(|)  cosX  —  cos<j)  sinX  —  sin(j) 


This  matrix,  an  array  of  direction  cosines,  is  also  available  with  high  precision  in  the  LN93  output. 
3.2  Two  Approaches  to  Vector  Gravimetry 

There  are  now  two  approaches  to  estimating  the  gravity  disturbance  vector.  The  first  approach 
uses  an  independent  measurement  of  position  (or  velocity),  such  as  from  GPS,  to  compare  against 
the  INS-indicated  position  (or  velocity).  The  difference  is  due  to  many  sources  of  error,  including 
the  fact  that  only  normal  gravity,  not  the  actual  gravity  vector,  was  used  in  the  integration  of  (7). 
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The  total  error  in  indicated  velocity  satisfies  a  linear  differential  equation  that  is  readily  derived 
(using  a  linearity  assumption  for  the  error)  from  (4),  as  follows: 

is**  =  -  5«  +  O  v«  -  «  +  O  5v”  +  aa-  +  r”  Sp"  +  8g»  (9) 

where  5p'*'  is  the  gravity  gradient  tensor  multiplied  by  position  errors  in  the  w-frame  (the 
horizontal  components  of  the  resulting  vector  are  usually  negligible).  The  acceleration  error,  6a , 
further  consists  of  sensor  errors,  5a*’ ,  and  a  Coriolis  term  due  to  orientation  errors,  t|f'^ ,  resulting 
from  the  erroneously  computed  transformation  matrix,  : 

5a’*'  =  C^5a*’  +  a’^x\|r’*'  (10) 

These  orientation  errors  satisfy  their  own  differential  equation  given  by  (see  Britting,  1971) 

¥“’  =  -<xv”'-CJ'5<oi+8«)|“„  (11) 

where5<»yj  are  the  gyro  errors.  We  note  that  8a>j^  and  5a)j^  in  (9)  and  (1 1)  are  also  functions  of 
errors  in  the  velocity.  Thus,  together,  (9)  and  (11)  yield  a  system  of  linear  differential  equations 
for  velocity  and  orientation  errors.  They  describe  a  model  for  the  dynamics  of  these  errors. 

The  sensor  errors  and  the  gravity  errors  each  have  their  own  dynamics.  The  accelerometer 
error  may  comprise  a  bias  (random  constant),  a  scale  factor  error  (random  constant),  correlated 
noise  (e.g.,  a  first-order  Gauss-Markov  process),  and  white  noise.  Similarly,  the  gyro  error  may 
consist  of  a  drift  (random  constant),  a  scale  factor  error  (random  constant),  correlated  noise  (first- 
order  Gauss-Markov  process),  and  white  noise.  Finally,  the  gravity  error  should  be  modeled  as 
some  kind  of  correlated  process.  Usually,  a  second-  or  third-order  Gauss-Markov  model  is  used. 
Including  these  dynamics  into  the  basic  dynamics  system  for  the  INS  errors  augments  the  system 
of  linear  differential  equations  significantly.  The  details  of  this  are  not  given  here;  they  may  be 
found  (for  the  case  of  n-frame  coordinatization)  in  (Jekeli,  1995)  and  various  other  articles 
referenced  therein. 

The  basic  idea  is  to  use  a  Kalman  filter  to  distribute  the  position  or  velocity  information, 
coming  from  an  external  source,  such  as  GPS,  in  an  optimal  way  so  as  to  estimate  the  various 
errors  of  the  system.  One  of  these  will  be  an  estimate  of  the  gravity  error,  that  is  the  gravity 
disturbance  vector.  Again,  the  details  of  the  formulations  of  the  Kalman  filter  and  all  associated 
models  are  not  pursued  here.  Only  a  schematic  diagram  is  shown  in  Figure  6  to  illustrate  the 
essential  processing  technique.  All  details  related  to  the  processing  of  GPS  phase  data  to  obtain 
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positions  (or  velocities)  are  also  omitted. 


apply  corrections 


double-difference 
phase  processing 


Figure  6:  Traditional  approach  to  estimate  gravity  disturbance  from  INS/GPS. 

We  note  that  the  estimation  of  the  gravity  disturbance  vector  with  this  method  absolutely 
requires  that  it  be  assigned  a  stochastic  model.  This  can  be  a  serious  drawback  to  the  method  since 
the  result  may  depend  critically  on  the  type  of  model  chosen  (see  Jekeli,  1995).  Furthermore,  only 
the  horizontal  components  of  the  gravity  disturbance  vector  can  be  estimated  this  way  since  the 
integration  of  the  inertial  accelerations  is  quite  unstable  in  the  vertical  (this  is  a  well-known 
phenomenon  that  precludes  free-inertial  navigation  in  the  vertical  for  periods  longer  than  a  few 
minutes). 

The  second  approach  is  fundamentally  different  from  the  first.  In  this  case  no  model  for  the 
gravity  disturbance  is  required,  at  least  in  theory.  Instead  of  time-integrating  the  IMU  data  to 
obtain  positions  (or  velocities),  the  GPS  positions  are  time-differentiated  to  obtain  accelerations. 
These  are  then  compared  to  the  accelerations  coming  from  the  IMU’s.  One  still  has  to  account  for 
Coriolis  terms  due  to  rotating  frames,  and  the  estimation  of  IMU  errors  can  still  be  done  with  a 
Kalman  filter,  but  the  gravity  disturbance  is  the  result  simply  of  a  final  differencing  between  two 
type  of  accelerations,  fundamentally  as  alluded  to  by  equation  (1).  This  method  is,  in  fact, 
precisely  the  way  airborne  scalar  gravimetry  is  performed.  Gravimeter  measurements  are  corrected 
for  the  Eotvos  effect  and  differenced  with  vertical  accelerations  determined  from  GPS,  or  a  laser 
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altimeter,  or  the  like. 

In  the  case  of  the  LN93,  where  raw  acceleration  data  are  not  available,  it  is  still  possible  to 
derive  a  method  based  on  this  principle.  We  simply  take  the  indicated  velocities  and  transform 
them  back  into  accelerations  with  an  appropriate  numerical  differentiation  algorithm.  From  (1)  and 
(2),  we  have 


5g^  =  Cj^x‘-a'*'-Y^ 


(12) 


where  Cf'  is  the  transformation  matrix  from  the  i-frame  to  the  w-frame.  Substituting  (7)  then 
yields 


8g-=crx‘-[^v;;'+{(B;;+<)xv 


(13) 


With  =  Cg  Cj  and  (5)  we  rewrite  this  as  follows 


6g^  =  CCfi‘ 


(14) 


where  we  note  that  C*  x‘  x® .  Also,  the  transformation  matrix,  Cf ,  and  the  angular  rate,  , 
can  be  computed  without  error  from  Earth’s  rotation  rate. 

It  is  possible  to  determine,  independently,  the  total  acceleration  Cj  x  from  GPS  (the  details, 
again,  are  omitted  as  being  outside  the  scope  of  the  present  report).  The  transformation  matrix, 

,  is  provided  as  output  by  the  LN93  computer  with  high  precision  (although  it  also  contains  the 
same  quantization  error  demonstrated  in  Figure  4,  this  will  affect  the  acceleration  computation  at 
the  part  in  lO”^  level  -  about  0.1  mgal).  The  angular  rate,  ,  is  given  by  (6)  (a  more  accurate 
form  that  accounts  for  ellipsoidal  effects  can  be  derived),  and  the  w-frame  velocities  are  provided 
by  the  LN93  with  high  precision.  It  remains,  then,  to  differentiate  these  velocities  numerically, 
and  all  quantities  on  the  right  and  side  of  (14)  are  thus  obtained.  Therefore,  it  is  possible  to 
determine  the  gravity  disturbance  from  indicated  LN93  velocities  and  direction  cosine  matrix  and 
from  GPS  accelerations.  The  errors  arising  form  the  IMU’s  must  be  treated  in  some  way.  This 
could  be  done  in  an  optimal  manner  as  discussed  below  for  the  case  of  raw  data  processing. 

In  the  case  that  raw  accelerometer  and  gyro  data  are  available,  as  from  CfM’s  LNIOO,  one  also 
has  two  data  processing  options  available  to  estimate  the  gravity  disturbance  components.  The 
traditional  approach  will  not  be  described  since  it  is  identical  to  the  case  above,  except  that 
positions  can  now  be  determined  independently  of  the  navigation  computer,  thus  ensuring 
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maximum  computational  precision  (i.e.,  no  artiHcial  quantization  effects).  Only  the  second  method 
based  on  a  comparison  of  accelerations  will  be  discussed  in  some  detail.  In  particular,  the 
processing  of  raw  gyro  and  accelerometer  data  will  be  described  in  Section  3.3  following  the 
essential  formula  for  gravity  disturbance  computation. 

This  second  approach  to  gravity  disturbance  vector  estimation  in  the  general  sense  has  been 
considered  by  Jekeli  and  Garcia  (1997),  as  well  as  by  Schwarz  and  Li  (1996).  And,  as  already 
mentioned  it  is  completely  analogous  to  traditional  scalar  gravimetry,  but  extended  to  three 
dimensions.  The  essential  idea  is  to  combine  equation  (2)  and  equation  (1)  as  follows: 


8g"  =  g"-r' 


(15) 


where  Cj,  is  the  transformation  to  the  i-frame  from  the  body  frame  (b-frame),  in  which  the  inertial 
accelerations,  ,  are  measured.  The  transformation  is  obtained  directly  from  the  gyros  of  the 
INS,  as  shown  below.  This  “direct  differencing”  method  is  much  simpler  computationally  than  the 
process  indicated  in  Figure  6,  and  it  uses  data  directly  from  the  gyros  and  accelerometers,  rather 
than  integrating  these  data  to  the  level  of  positions.  On  the  other  hand,  the  GPS  position  data  must 
now  be  elevated  to  the  level  of  acceleration.  This  means  that  noise  in  GPS  position  data  is 
amplified  at  high  frequencies  and  requires  careful  filtering. 


3.3  Algorithms  to  Process  Raw  Gyro  and  Accelerometer  Data 

The  IMU  produces  accelerometer  and  gyro  pulses:  5vj  and  BQf ,  respectively;  where  5v|  is  a 
vector  of  increments  in  the  sensor-frame  velocity  generated  by  the  three  accelerometers,  and  50|  is 
a  vector  of  increments  in  the  sensor-frame  angles  generated  by  the  three  gyros.  With 
corresponding  time  increments,  8t ,  we  have 


8v|  = 


aWdt  , 


& 


80,= 


fi»is(t)  dt 


tt 


(16a,b) 


where  a®  is  the  acceleration  vector  in  the  s-frame  (sensor  frame)  and  is  the  angular  rate  with 
respect  to  the  i-frame  (inertial  frame),  also  in  the  s-frame.  Note  that  the  integration  takes  place  in 
the  s-frame  which  is  rotating  with  respect  to  the  inertial  frame.  Also,  we  do  not  assume  that  the 
sensor  and  body  frames  are  equivalent.  In  fact,  in  determining  accelerations  we  must  distinguish 
between  the  locations  of  the  INS  and  the  GPS  antenna.  This  difference  causes  the  so-called  lever- 
arm  effect.  The  body  frame  is  usually  associated  with  the  aircraft  and  we  will  assume  the  GPS 
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antenna  to  be  the  origin  of  that  frame.  The  s-frame  has  its  origin  located  at  some  fiducial  point  on 
the  INS,  as  specified  by  the  manufacturer.  The  lever-arm  effect  is  discussed  later  in  Section  3.3.3. 

To  determine  the  acceleration  in  the  i-frame,  it  is  necessary  to  solve  for  the  acceleration  in  the 
sensor  frame  using  the  5v(  ’s,  and  then  to  rotate  this  into  the  i-frame  using  the  orientation 
information  contained  in  the  50( ’s.  The  rotation  matrix,  Cg ,  that  transforms  a®  into  a* , 
according  to  a®  =  a®  is  determined  by  integrating  a  differential  equation  relating  the  angular  rates 
to  the  rotation  matrix: 


(17) 


where  f2®s  is  a  skew-symmetric  matrix  with  components  of  (»?,  as  its  off-diagonal  elements: 


0  -©3  ©2 

©3  0  -©1 

-©2  ©1  0 


(18) 


where  =  [©j,  ©2  ©3]^ . 

It  is  far  easier  to  formulate  an  equivalent  differential  equation  in  terms  of  quaternions.  Let 
q  =  (a,  b,  c,  d)’^  be  a  four-vector  of  time-dependent  quantities  such  that  q  satisfies  the  following 
linear  system  of  differential  equations 

q  =  -jAq  (19) 


where  A  is  a  4x4  skew-symmetric  matrix  of  time-dependent  angular  rates: 


0  ©1  ©2  ©3 

-©2  0  ©3  -©2 

—©2  0 

—©3  0)2  — ©2  11 


It  can  be  shown  that 


a^  +  b^-c^-d^  2  (be -tad)  2(bd-ac) 

2  (be -ad)  a^-b^-tc^-d^  2(cd-t-ab) 

2  (bd  +  ac)  2  (cd  -  ab)  a^  -  b^  -  c^  -t-  d^ 


(20) 


(21) 
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The  elements  of  q  are  known  as  quaternions.  The  integration  of  q  is  done  incrementally  over  a 
time  interval,  such  as  5t ,  or  a  multiple  of  this;  and,  this  gives  the  elements  of  C[  for  each  step. 

3.3. 1  A  Second-Order  Algorithm 

Since  the  time  increment  for  an  IMU  is  usually  very  small;  e.g.,  the  data  rate  may  be  256  Hz, 
so  that  5t  =  0.00390625  s ,  it  may  be  considered  that  is  constant  over  this  time  interval.  It  is 
not  necessary  to  assume  this;  but  then  the  integration  formulas  become  more  complicated,  as  seen 
in  the  subsequent  section.  For  constant  angular  rates  over  5t ,  we  have  the  estimate,  from  (16b), 
for  the  step: 


S0|  =  ®L(tj)8t  (22) 

and,  with  matrix  A  assumed  constant  over  this  interval,  denoted  by  A| ,  an  analytic  form  of  the 
solution  to  (19)  can  readily  be  found.  It  is 

q(t)  =  ei^'^‘"**-'^q(t|_i)  ,  0<t-tj_,<5t  (23) 

which,  in  series  form,  yields  the  following  iterative  solution: 

q(t, )  =  (l  i  A*  St  +  i  At  St^  +  ^  A?  8t^  + ...)  i(t|_  j) ,  f  =  1,2,...  (24) 

where  a  single  initial  condition,  q(to) ,  must  be  specified.  A  corresponding  closed  form  may  be 
found  by  noting  that  the  powers  of  the  matrix  Aj  are  either  diagonal  (even  powers)  or 

A 

proportional  to  Aj  itself  (odd  powers).  One  has: 


q(t|)  = 


cos(^  1 80, 1)  I  +  sin(^  1 80,  |)  B, 


qdt  _  i) 


where  1 80,|  =  ^  80^80,  and 


(25) 


B,  =  A,8t  (26) 

The  solution  (24)  is  given,  in  principle,  with  two  errors,  an  algorithm  error  (the  solution  (23) 
assumed  constant  A-matrix)  and  a  data  error,  both  due  to  the  assumption  of  constant  angular  rates. 
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However,  the  data  error,  in  a  sense,  is  fictitious  since  the  data  are  not  the  angular  rates,  as 

^  _ 

approximated  by  A{ ,  but  the  integrals  given  by  (16b).  Thus,  let 


B*  = 


A(t)  dt 


(27) 


Using  Bj  instead  of  Aj  5t  in  (24)  yields  the  actual  algorithm 
5(t, )  =  (l  +  i  Bf  + 1  i  B?  + ...)  5(t,_  i) .  C  =  1,2,... 


(28) 


This  formula,  however,  still  contains  algorithm  error. 

This  algorithm  error  is  of  the  order  5t^ .  Indeed,  let  the  true  quaternion  be  represented  over  the 
integration  interval  by  the  series 


q(tj)  =  q(t(_  i)  +  q(t{_  i)  5t  +  q(t{_  i)  5t^  +  ^  q(t(_  j)  8t^  + 


(29) 


Substituting  (19)  yields 


q*  -  1  St + i(At_  1  +  8t^ 

+  ^(A{_i  +  A£_i  A|-i  +5A|_i  A(_i  +iAj_ij8t^  +  ...]  qj_i 


(30) 


where  all  quantities  are  true,  not  approximate,  and  subscripts  denote  the  time  epoch  at  which  they 
are  defined.  Similarly,  from  (27)  we  may  write  in  terms  of  true  quantities: 


Bj  =  A|_i  5t  +  ^  A(_|  5t'^  +  ^  A{_i  8t^  +. 


2  .  1  X 


(31) 


Substituting  (31)  into  (28)  and  comparing  with  (30),  it  is  easily  verified  that  the  algorithm  error  is 


it  -  q{  =  [i(A{- 1  A(_  1  -  A(_  1  A{_  i)  5t^  +  ...]  qf_  j 


(32) 


assuming  the  previous  error,  qj-i  -  q(_  i ,  yields  higher-order  terms. 
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It  can  be  shown  that  the  first  term  in  (32)  vanishes  if  the  angular  rate  vector  does  not  change 
direction  during  the  integration  interval;  i.e.,  in  that  case  the  A-matrix  and  its  first  derivative 
commute.  That  is  why  the  algorithm  (28)  is  said  to  have  commutativity  error,  or  coning  error. 
In  other  words,  if  the  angular  rate  vector  does  change  direction,  the  system  is  “coning”. 

3.3.2  A  Third-Order  Algorithm 

One  can  design  higher-order  algorithms  in  a  number  of  ways.  One  is  to  do  a  numerical 
integration  of  (19)  based  on  a  model  for  the  angular  rates.  The  Runge-Kutta  (R.K.)  algorithm  is 
one  such  integrator;  and,  in  turns  out  that  the  third  order  R.K.  method  with  a  linear  model  for  the 
angular  rates  yields  an  algorithm  for  the  quaternions  with  fourth-order  algorithm  error.  The  third- 
order  R.K.  method  requires  that  the  function  being  integrated  is  evaluated  at  either  end  of  the 
integration  interval  and  half-way  in  between.  Therefore,  in  this  case  the  integration  interval  is 
twice  the  data  interval: 

At  =  2  5t  (33) 


To  simplify  the  notation,  let  ca?g(t)  s  co(t) ,  and  assume  that  over  the  integration  interval: 

®(t)  =  0^-2  +  ^-2(t-te-2)  +  O(At2)  ;  |t-t,_2|<At  (34) 


where  the  subscripts  denote  evaluation  of  the  (Uiae)  quantity  at  the  corresponding  epochs  spanning 
intervals  5t ,  not  At .  Substituting  (34)  into  (16b)  yields 


60|_ 


V 

CB(t')  df  =  2  5t  +  i  Q3^_ 2  6t2  +  0(At3) 

i-2 


(35) 


50. 


(0(0  df  =  (aj_2  8t  +  J  ^_2  +  O(At^) 


(36) 


Solving  for  (iin_2  and  ^_2 1  it  is  readily  verified  that 
**^-2  =  ^(3  50{_  1  -  50j)  +  O(At^) 


(37) 
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oq^_2  =  (80f  —  80{_  j)  +  0(At)  (38) 

Using  these  expressions  in  (34)  yields  for  three  consecutive  epochs: 

QH_2At  =  380{_i-5ej  +  O(At^)  (39) 

OQj.j  At  =  80t_i  +  80{  +  O(At^)  (40) 

OH  At  =  3  80J  -  80, _  1  +  0(At3)  (41) 

Therefore,  we  have  the  following  observed  quantities,  which  are  accurate  in  terms  of  the  model  to 
second  order: 

an_2At  =  3  80{_i-80,  (42) 

ca,_i  At  =  80,_i +  80,  (43) 

^  At  =  3  80,-80,_i  (44) 

With  the  Taylor  expansions  (35)  and  (36)  substituted  into  the  right  sides,  we  also  obtain 
0fl^_2  At  =  cin_2  At  +  0(At^) 

c^_jAt  =  tH^_2At  +  -jOH_2At^  +  0{At^)  (46) 

on  At  =  ai,_2  At  +  ^_2  At^  +  0(At^)  (47) 

The  third-order  R.K.  algorithm,  for  q  =  f(q,t) ,  is  given  by 

i  =  i-  2  + 1  (^qo  +  4  Aq  1  +  Aq2)  (48) 

where 
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Aqo  =  f(q{-2.  tt_2)At 

(49) 

Aqi  =  f(q(_2  +  |Aqo,tj_,jAt 

(50) 

Aq2  =  ^qt-2  +  2  Aqi  -  Aqo,  t{)  At 

(51) 

This  iterative  integration  algorithm  provides  an  approximate  solution  at  intervals  of  At ;  but  one 
could  do  a  second  iteration,  using  the  same  data,  on  the  in-between  points  -  it  is  just  a  matter  of 
which  starting  value  is  chosen.  Note  that  only  one  starting  quaternion  is  required;  the  R.K. 
algorithm  is  a  single-step  method. 

From  (19),  we  have  f(q,  t)  =  A(t)  q/2  ;  and  using  the  observed  values  (42)  -  (44)  for 
A(t|_2),  A(t|_  i),  A(t|) ,  equations  (49)  -  (51)  become  with  Bj  =  Aj  At  (from  (42)  -  (44)): 

Aqo  =  ^B{_2qj_2 

(52) 

Aqi  =  ^B(_i  |l-i-iBj_2jq{_2 

(53) 

Aq2  =  ^  B{  +  Bj_  1  ^1  1  B,_2j  -  ^  B,_2  qf.2 

(54) 

Substituting  these  into  (48)  and  simplifying  the  result  yields 


qt= 


1+ 


i  (b,  +  4  B,. .  +  B,_  j)  +  ±  (l  +  ^  §,)  B,_ ,  B,_  J  +  ±  B,  (b,.  ,  -  i  B..^]]  i. 


(55) 


By  substituting  (45)  -  (47)  into  (55),  it  is  readily  verified,  with  a  comparison  to  (30),  that 

q{-qe  =  0(5t^)qj_i 

That  is,  equation  (55)  is  a  third-order  algorithm;  the  algorithm  error  is  of  fourth  order. 

Higher  order  algorithms  can  be  derived,  but  they  will  require  larger  integration  intervals  per 
step  (if  the  basic  data  interval,  5t ,  remains  the  same). 

The  initial  value  for  the  iteration,  qo  (or,  q, ),  is  obtained  from  the  initialization  (alignment)  of 
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the  IMU,  typically  provided  by  the  navigation  computer.  Or  one  can  determine  the  initial  attitude 
externally  using  (zero)  velocity  updates  in  a  Kalman  filter  using  standard  procedures.  If  the  result 
is  an  initial  rotation  matrix  the  following  inverse  relationship  to  (21)  can  be  used  to  get  the  initial 
quaternion: 


^  “  2  Id  (^s)2,2 

^  ^  ({^^2,3  “  (^03,2) 

(57) 

^~4a  ((^s)3,i  ~ 

^  ((^s)  1.2  “  (^5)2,1) 


These  equations  are  readily  verified  since  C*  is  an  orthogonal  matrix,  where 

Ci  =  Ci,C^  (58) 

and  n  denotes  the  navigation  frame  (n-frame).  C"  is  obtained  during  the  initialization  phase  and 
C},  can  be  determined  from  the  known  position  of  the  system  at  the  initial  time. 

Finally,  it  is  noted  that  the  quaternions  must  satisfy  the  following  property: 


a2  +  b^  +  c2  +  d^=  1  (59) 

On  the  other  hand,  the  computed  quaternions,  q ,  from  any  algorithm  may  not  satisfy  this 
constraint.  This  is  due  to  numerical  round-off  error  as  well  as  model  error  and  will  cause  the 
computed  transformation  matrix  to  become  non-orthogonal.  To  avoid  this  it  is  advisable 
(practically  necessary)  to  “re-orthogonalize”  the  computed  quaternions  at  each  iteration  according  to 

9  ^-^5  (60) 

vq  q 


^  ^  a 

which  will  ensure  that  q  q  =  1  and  that  is  orthogonal. 
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3.3.3  Lever-Arm  Effect 

Consider  the  sensor  frame  to  be  attached  to  the  IMU,  as  before.  In  order  to  compare  the 
inertial  accelerations  determined  at  the  IMU  with  the  inertial  accelerations  determined  by  GPS  at  the 
GPS  antenna,  it  is  necessary  to  correct  for  the  lever-arm  effects  of  the  IMU  with  respect  to  the 
antenna  (or  vice  versa).  Without  loss  in  generality,  one  may  assume  that  the  b-frame  and  s-frame 
are  paraUel;  but  their  origins  are  offset  by  the  lever  arm,  b** ,  the  vector  of  coordinates  of  the  sensor 
(IMU)  in  the  b-frame  (Figure  7). 

We  have  the  following  relationship 

ZS  V*  —  V* 

accel  *body  (61) 

where  all  coordinates  are  in  the  i-frame.  The  acceleration  of  the  b-frame  can  be  expressed  in  terms 
of  the  sensed  accelerations  in  the  s-frame  by  applying  the  the  Coriolis  law  to  (61),  noting  that 
b  =  0  (the  two  frames  are  rigidly  connected,  by  assumption): 

*accel~*body'*‘'^  =*body‘*"Ci,0)|jjXb*’ 


Differentiating  again,  we  obtain 


*accel  =  Xbody  +  (^b  tojj,  +  Cj,  ©y  X  b**  =  x{^y  +  Cj,  ©J,  X  b*’  +  C},  ©J,  X  (©jj^  X  b**) 
Coordinatizing  this  in  the  b-frame  using  C?* ,  we  have  with  x=a'  +  g*  from  (1): 

»accel  =  ^accel  =  +  E^y  “  B^cel  +  ®ib  X  b**  +  ©Jj,  X  (©!’j^  X  b^) 


Where  it  is  assumed  that  C*  =1.  This  says  that  the  acceleration  sensed  in  the  s-frame,  but 
coordmatized  in  the  b-frame  to  which  it  is  rigidly  attached,  is  equal  to  the  acceleration  of  the  b- 
frame  plus  various  reaction  forces  supplied  by  the  rigid  support  of  the  s-frame  to  the  b-frame.  The 
first  is  due  to  the  gravitational  difference  between  locations  of  the  accelerometer  and  case  frame 
origins,  and  the  others  are  associated  with  the  rotation  of  the  case  in  the  i-frame.  g**  is  the 
gravitation  vector  at  the  indicated  point  in  the  b-frame;  the  only  component  of  concern  is  the  vertical 
component  which  varies  by  about  -  0.3086  mgal/m . 
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Alternative  formulas  for  the  lever-arm  effect  from  a  practical  viewpoint  can  be  derived  for  the 
case  of  twice  differentiating  lever-arm  positions  in  the  i-frame: 

a^cel  =  +  gtody  ~  Saccel  (65) 

and  for  the  case  of  once  differentiating  lever-arm  velocities,  also  in  the  i-frame: 

®accel  “  ®body  Sbody  ~  Sacccl"^  )  (66) 

It  is  readily  verified  that  (65)  and  (66)  are  equivalent  to  (64).  For  the  practical  computations,  in 
both  cases,  it  is  necessary  to  know  the  transformation,  C{, ,  from  the  body  frame  to  the  inertial 
frame.  For  the  latter  case,  also  the  angular  rate  of  rotation,  co-J, ,  must  be  known.  These  quantities 
are  determined  from  the  gyro  data  using  an  algorithm  such  as  described  in  Section  3.3.2.  Note  that 
the  transformation  of  the  lever-arm  to  the  i-frame  must  be  done  before  implementing  a  time- 
differentiation  algorithm.  Experiments  have  shown  that  these  formulas  are  more  stable  and 
accurate  than  (64)  since  the  latter  requires  a  numerical  time-differentiation  of  the  angular  rate,  . 
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3.3.4  Final  Computation 

To  compute  the  sensed  accelerations  form  the  accelerometer  output  data  (16a),  we  may  use  a 
first-order  estimate: 

at=^(5v,^.,+5v{)  (67) 

Having  determined  the  quaternions,  q(t|) ,  using  (55),  one  computes  C‘(tj)  from  (21).  At  the 
same  time,  the  lever-arm  effects  are  applied  to  the  measured  accelerations  (67)  according  to  (65),  or 
(66).  The  normal  gradient  of  gravitation  should  be  used  if  there  is  a  significant  vertical  component 
in  the  lever-arm.  The  result  is  the  b-frame  accelerations,  a‘’(t{) ,  at  the  GPS  antenna.  Then,  since 
Cb(tj)  =C*(t|) ,  there  is 

a‘(ti)  =  C|,(t|)  a*’(tt)  (6g) 

where  a'(t()  is  the  inertial  acceleration  in  the  i-frame.  This  can  be  compared  directly  to  the  GPS- 
derived  kinematic  acceleration  in  the  i-frame,  x* ,  to  get  the  gravitation  vector  in  the  i-frame: 

g‘  =  x‘-a‘  (69) 

It  remains  only  to  transform  the  gravitation  vector  into  more  appropriate  coordinates,  the  n- 
fiame: 


8"  =  C?gi  (70) 

where,  again,  C|j  can  be  determined  from  the  position  of  the  GPS  antenna.  The  accuracy  of  Cjj 
need  be  only  1  part  in  10^  to  get  mgal  accuracy  in  the  gravitation  vector. 

This  scheme  so  far  does  not  account  for  the  IMU  errors.  One  method  to  estimate  these  errors 
is  to  use  a  Kalman  filter  on  a  model  for  the  accelerations,  where  the  GPS-derived  accelerations 
serve  as  updates  to  the  filter.  The  residual  between  the  observed  GPS  accelerations  and  the 
adjusted  GPS  accelerations  would  indicate  the  gravity  disturbance.  In  this  case  no  stochastic 
model  for  the  gravity  disturbance  is  required  as  in  the  case  of  Kalman-filtering  INS  positions.  A 
schematic  of  this  procedure  is  given  in  Figure  8. 
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Figure  8:  Alternative  Data  Processing  for  Vector  Gravimetry 


4  Tests  and  Analysis  with  the  LN93  and  LN100 

Much  of  the  difficulties  surrounding  the  use  of  the  LN93  centered  on  understanding  the  user 
interface  software  supplied  with  the  system.  The  software,  written  in  Pascal,  is  authored  by  Litton 
and  is  called  the  User-Friendly  Bus  Control  (UFBC)  software.  It  is  written  for  the  desk-top 
personal  computer  and  a  Pentium  based  system  is  recommended.  Furthermore,  since  a  high 
volume  of  data  from  the  INS  must  be  written  to  data  storage  veiy  quickly,  a  RAMDRTVE  must  be 
set  up  in  RAM;  at  least  10  MB  of  memory  must  be  available  for  the  RAMDRTVE.  Some  details  on 
the  individual  components  of  the  interface  are  documented  in  the  report  by  Humphrey  and 
Kawakami  (1996). 

Figure  9  shows  a  schematic  of  the  system  interrupts  that  are  used  to  control  the  timing  and  data 
requests  of  the  1553  bus  controller.  The  timing  of  the  INS  data  is  crucial  since  the  difference 
between  GPS  and  INS  clocks  can  amount  to  significant  lags.  The  entire  data  retrieval  from  the 
INS  is  governed  by  GPS  time.  A  signal,  called  the  PPS  (pulse-per-second),  is  sent  by  the  GPS 
receiver  through  the  timer  card  to  the  PC.  The  PPS  is  synchronized  to  GPS  time  to  within  a 
millisecond,  or  better,  depending  on  the  receiver.  Receipt  of  the  PPS  by  the  PC  initiates  a 
sequence  of  events.  First,  a  “no-data”  message  is  sent  to  the  INS,  causing  a  bit-failure  in  the 
system  indicating  a  timer  problem.  This  causes  the  INS  clock  to  be  reset  to  zero.  This  clock  is  a 
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counter  that  fills  a  16-bit  register  with  64  microsecond  counts.  After  about  4.2  s,  the  counter 
“turns  over”  to  zero;  thus  a  reset  from  the  GPS  PPS  will  yield  an  unambiguous  time  tag  from  the 
INS  between  zero  and  one  second  for  each  retrieved  data  item.  It  is  only  necessary  to  identify  the 
PPS  in  terms  of  absolute  time  from  the  GPS  receiver. 


Figure  9:  Interrupts  that  control  the  timing  of  INS  data  requests. 

The  next  event  triggered  by  the  receipt  of  the  PPS  is  the  start  of  a  clock  on  the  timer  board  that 
counts  off  0.05  seconds,  where  each  count  serves  to  identify  the  next  request  of  data  from  the  INS 
through  the  1553  bus.  It  is  noted  that  the  time  tags  of  the  INS  data  refer  to  the  time  of  validity  of 
the  data  computed  by  the  INS,  for  example,  the  velocity  or  position  coordinates.  These  time 
epochs  reflect  the  time  at  which  the  raw  data  were  sampled  to  compute  the  position  or  velocity. 
Since  there  is  a  signiHcant  time  (typically  <  0.1  s)  associated  with  these  computations,  the  time  tag 
at  the  epoch  of  the  PPS,  and  for  two  epochs  beyond  this,  actually  refers  to  the  previous  PPS.  This 
is  important  to  note  for  the  proper  interpretation  of  the  time  tags  of  the  INS  data  relative  to  the  GPS 
time. 
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It  is  also  noted  that  the  bit-failure  caused  by  the  no-data  message  prevents  the  INS  from 
optimally  initializing  the  INS  and  entering  into  the  NAV-READY  mode.  To  circumvent  the  bit- 
failure  on  initialization,  a  different  bus  controller  file  (that  does  not  send  this  message)  is  used  to 
allow  the  INS  to  achieve  the  NAV-READY  mode.  Once  NAV-READY  is  achieved,  a  new  bus 
controller  file  must  be  initiated  by  the  user. 

The  UFBC  software  was  modified  to  allow  the  option  of  running  the  INS  without  GPS  time 
synchronization.  This  proved  useful  in  tests  of  the  UFBC  software  and  in  laboratory  tests  where 
only  the  navigation  error  characteristics  of  the  LN93  were  to  be  investigated. 

)  Some  difficulties  that  plagued  the  LN93  system  integration  included  the  failure  to  adapt  the 

UFBC  software  to  retrieve  and  save  also  the  GPS  data  from  the  receiver.  This  may  be  due  to  the 
high  I/O  rates  that  are  imposed  on  the  system,  causing  failures  in  saving  any  of  the  required  data. 
It  was  judged  that  the  GPS  receiver  should  save  the  GPS  data  independently  of  the  INS,  where,  of 
course,  the  same  PPS  time  tags  are  applied  to  the  GPS  data.  Other  problems  that  needed  to  be 
overcome  for  practical  system  integration  concerned  providing  appropriate  power  from  12-VDC 
batteries,  or  from  28-VDC  aircraft  power,  to  the  28  VDC  input  power  of  the  INS,  the  120V AC 
input  power  of  the  PC,  and  the  12VDC  input  power  of  the  GPS  receiver. 

Initial  tests  of  the  LN93  were  done  in  the  laboratory  in  order  to  understand  the  INS  output  data 
and  check  out  the  general  operation  of  the  system.  These  test  were  followed  by  road  tests,  where 
the  INS  was  strapped  in  the  back  of  a  van  and  driven  over  a  period  of  about  45  minutes  in  the  area 
of  Columbus.  Figure  10  shows  the  closed  course  followed  for  three  such  tests.  One  course  could 
not  be  completed  due  to  a  power  failure  just  prior  to  returning  to  the  start  point.  The  misclosure, 
calculated  as  the  total  distance  between  the  true  and  indicated  end  points  for  each  course,  on 
average,  was  670  m.  This  is  well  within  the  0.8  nmi/hr  =1.4  km/hr  speciOed  for  the  system;  and 
it  shows  that  the  LN93  is  performing  as  expected,  or  better. 
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During  these  road  tests,  GPS  data  were  also  collected  using  a  Trimble  4100  SSI  receiver,  but 
the  results  were  rather  poor  due  to  frequent  signal  outages  caused  by  many  overpasses  and  other 
obstructions.  It  was  not  deemed  worthwhile  to  attempt  further  tests  with  the  integrated  GPS/INS 
on  the  van  at  this  time  because  of  the  difficulties  of  GPS  data  collection,  the  frequent  power 
problems  that  were  encountered,  and  the  difficulties  of  adapting  the  UFBC  software  to  the  specific 
applications  at  hand. 

Finally,  some  effort  was  undertaken  to  conduct  an  airborne  test  of  the  LN93.  In  collaboration 
with  the  Ohio  Department  of  Transportation  (ODOT),  a  cradle  was  constructed  for  the  INS  that 
could  be  slipped  onto  the  seat  rails  of  a  small  twin-engine  photogrammetric  aircraft.  One  test  was 
conducted  but  ended  in  failure  because  of  a  power  interruption  to  the  system  during  landing  and 
UFBC  data  handling  problems.  Future  tests  may  be  resumed  on  an  opportunity  basis. 

The  Center  for  Mapping  (CfM)  at  OSU  has  developed  an  Airborne  Integrated  Mapping  System 
(AIMS)  that  includes  a  Trimble  GPS  receiver  and  a  Litton  LNIOO  INS  (Grejner-Brzezinska  and 
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Phuyal,  1998).  This  system  has  undergone  numerous  flight  tests  and  the  data  from  one  of  the 
early  tests  were  analyzed  to  determine  the  viability  of  this  system  for  vector  gravimetry.  Ground 
tests  conducted  by  CfM  similar  to  those  described  above  indicated  that  also  the  LN 100  was 
performing  at  the  specified  level  of  0.6  nmi/hr  =  1  km/hr. 

The  test  flight  that  was  analyzed  was  conducted  in  the  St.  Louis,  Missouri,  area;  and  the 
indicated  (i.e.,  free-inertial)  INS  and  GPS  trajectories  are  shown  in  Figure  11.  The  misclosure  of 
the  INS  trajectory  is  about  30  km  after  less  than  one  hour  of  flight,  probably  due  to  excessive 
errors  in  the  gyro  data  (clearly,  the  system  did  not  perform  within  specification  during  this  test). 
Grejner-Brzezinska  and  Phuyal  (1998)  report  that  for  this  test  flight  the  attitude  error  estimated 
using  a  Kalman  and  GPS  position  updates  had  standard  deviations  up  to  17  arcsec  (rms). 
Therefore,  this  flight  test  was  not  particularly  suitable  for  the  purpose  of  vector  gravimetry. 


Figure  11:  INS  and  GPS  trajectories  of  Flight  #3,  CfM,  3/6/97. 


However,  in  order  to  test  the  Kalman  filter  software  for  estimating  the  gravity  disturbance 
vector,  and  because  the  trajectory  included  a  relatively  straight  path,  the  data  were  analyzed.  The 
actual  deflection  of  the  vertical  along  the  flight  path  was  interpolated  from  a  2’x2’  DOV  data  base 
provided  by  the  National  Imaging  and  Mapping  Agency  (NIMA).  Figure  12  shows  the  relative 
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deflections  along  the  flight  path;  and  Figure  13  shows  the  estimated  deflections,  the  actual 
deflections,  and  their  differences.  The  Kalman  filter  that  was  used  to  estimate  the  deflections 
included  the  error  states  noted  in  Table  2,  as  well  as  a  correlated  noise  model  (first-order  Gauss- 
Markov)  for  the  accelerometers  with  variance  equal  to  25  mgal^  and  a  correlation  time  of  5  min. 
The  gravity  disturbance  model  was  a  third-order  Gauss  Markov  model  with  variance  equal  to 
480  mgal^  and  a  correlation  distance  of  about  26  km. 
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The  results  of  the  estimations  are  rather  poor  as  seen  in  Figure  13.  The  errors  in  the  estimated 
deflections  are  as  large  as  the  deflections  themselves.  This  is  very  likely  due  to  the  large  errors  in 
the  gyros.  However,  the  east  deflection  error  consists  primarily  of  a  trend  and  may  reflect  the  fact 
that  the  north  gyro  (whose  errors  contribute  to  the  east  deflection  error)  is  always  easier  to  calibrate 
during  INS  initialization  than  the  east  gyro.  We  also  see  from  Figure  1 1  that  most  of  the  position 
error  (at  least  initially)  is  in  the  north  direction,  again  pointing  to  large  east  gyro  errors.  Thus, 
while  the  numerical  results  of  this  test  are  not  satisfactory,  there  is  some  indication  that  we 
understand  the  reason  for  this,  which  in  itself  is  a  useful  conclusion. 


Figure  13.  a)  Estimates  of  deflection  components;  b)  true  deflection  components;  c) 
difference  between  estimates  and  true  deflections.  In  all  cases  the  solid  line  represents  the 
north  deflection  component.  Ordinate  values  in  units  of  arcseconds;  abscissa  in  seconds. 
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5  Summary  and  Conclusion 


This  report  describes  some  practical  experiences  with  the  LN93  and  LNIOO  inertial  navigation 
systems  in  the  area  of  airborne  vector  gravimetry.  No  conclusive  results  from  flight  test  data 
analyses  can  be  reported  this  time,  but  several  aspects  of  such  a  gravimetry  system  have  been 
elucidated.  The  two  inertial  navigation  systems  are  designated  by  the  manufacturer  as  medium-to- 
high  accuracy  systems,  indicating  that  the  measurements  are  of  sufficient  quality  that  when  < 

combined  with  geodetic  quality  GPS  receivers  they  should  be  capable  of  providing  the  precision 
needed  to  predict  the  gravity  vector  with  several  mgal  precision  (not  accuracy,  because  of  long¬ 
term  drifts  in  the  gyro  data).  The  two  systems  under  this  discussion  have  significantly  different 
data  output  capabilities.  In  each  case  special  data  processing  algorithms  must  be  developed  to 
make  optimum  use  of  the  data.  These  algorithms  are  described  in  detail  in  Section  3.  They 
differentiate  between  the  type  of  output,  whether  precise  velocities  and  direction  cosines  (LN93)  or 
raw  accelerometer  and  gyro  data  (LNIOO),  and  the  type  of  methodology  used  to  combine  the  INS 
and  GPS  data,  whether  by  time-integrating  the  INS  data  or  by  time-differentiating  the  GPS  data. 

Each  of  these  possible  processing  algorithms  carries  its  own  advantages  and  disadvantages. 

Using  the  velocities  determined  by  the  INS  navigation  computer  may  provide  the  best  available 
INS  velocities,  since  they  are  based  on  the  manufacturer’s  (in  this  case  Litton’s)  extensively  tested 
navigation  algorithm.  Using  the  raw  accelerometer  data  yields  more  flexibility  in  applying  different 
methodologies  to  combine  the  INS  and  GPS  data,  for  example,  allowing  them  to  be  more  tightly 
coupled.  Using  a  Kalman  filter  to  estimate  the  gravity  disturbance  components  from  time- 
integrated  IMU  data  is  founded  on  optimal  estimation  theory,  but  requires  a  stochastic  model  for 
the  gravity  disturbance.  Also,  the  vertical  component  cannot  be  estimated  since  the  IMU  data  are 
not  stable  when  integrated  in  the  vertical  dimension.  Using  a  straightforward  differencing  of  IMU 
and  time-differentiated  GPS  accelerations  to  obtain  the  gravity  disturbance  vector  is  a  natural 
extension  of  the  quite  successful  airborne  (scalar)  gravimetry  concept.  However,  a  unified 
approach  that  optimally  accounts  for  IMU  errors  has  yet  to  be  developed  and  tested. 

Several  system  tests  were  conducted  in  the  laboratory  and  in  the  field  to  determine  the  quality 
of  the  INS  and  to  begin  testing  the  algorithms  for  vector  gravimetry.  These  tests  showed  that  the 
INS’s  are  performing  as  expected.  However,  several  complications  in  physical  system  integration 
and  limited  availability  of  suitable  test  data  have  precluded  extensive  analysis  of  these  algorithms  at 
this  time.  It  is  anticipated  that  continued  development  of  data  processing  methodologies  and  more 
application-specific  airborne  testing  will  bear  fruit  in  these  areas.  The  simulation  studies  conducted 
by  numerous  authors  and  the  limited  test  results  shown  here  all  point  in  that  direction. 
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